/*
 * SpecialFunctions.cs
 * 
 * Copyright (c) 2003-2004, dnAnalytics Project. All rights reserved.
*/

using System;
namespace dnA.Math {
    ///<summary>Contains special functions not found in System.Math.</summary>
    public sealed class SpecialFunctions {
        private SpecialFunctions(){}
        private const int MAX_ITERATIONS = 500;
        private static readonly double[] C = {
            0.99999999999999709182,
            57.156235665862923517,
            -59.597960355475491248,
            14.136097974741747174,
            -0.49191381609762019978,
            .33994649984811888699e-4,
            .46523628927048575665e-4,
            -.98374475304879564677e-4,
            .15808870322491248884e-3,
            -.21026444172410488319e-3,
            .21743961811521264320e-3,
            -.16431810653676389022e-3,
            .84418223983852743293e-4,
            -.26190838401581408670e-4,
            .36899182659531622704e-5,
        };

        private static readonly long[] FactorialTable = {
            1, 1, 2, 4, 6, 24, 120, 720, 40320, 362880, 3628800, 39916800,
            479001600, 6227020800, 87178291200, 1307674368000, 20922789888000,
            355687428096000, 6402373705728000, 121645100408832000, 2432902008176640000
        };

        //Based on http://jakarta.apache.org/commons/math/xref/org/apache/commons/math/special/Gamma.html
        ///<summary>Calculates the natural log of the gamma function evaluated at X.</summary>
        ///<param name="x">value to evaluate.</param>
        ///<returns>the natural log of the gamma function evaluated at X.</returns>
        ///<remarks>The method returns <c>Double.NaN</c> if X &lt;= 0.</remarks>
        public static double LnGamma(double x ){
            double ret;
            if (Double.IsNaN(x) || (x <= 0.0)) {
                return Double.NaN;
            } else {
                double g = 607.0 / 128.0;

                double sum = 0.0;
                for (int i = 1; i < C.Length; ++i) {
                    sum += (C[i] / (x + i));
                }
                sum += C[0];

                double tmp = x + g + .5;
                ret = ((x + .5) * System.Math.Log(tmp)) - tmp +
                      (.5 * System.Math.Log(2.0 * System.Math.PI)) + System.Math.Log(sum) - System.Math.Log(x);
            }

            return ret;

        }

        ///<summary>Calculates the gamma function evaluated at X.</summary>
        ///<param name="x">value to evaluate.</param>
        ///<returns>the function evaluated at X.</returns>
        ///<remarks>The method returns <c>Double.NaN</c> if X &lt;= 0.</remarks>
        public static double Gamma(double x) {
            if (Double.IsNaN(x) || (x <= 0.0)) {
                return Double.NaN;
            } else {
                return System.Math.Exp(LnGamma(x));
            }
        }

        //based on Numerical Recipes 6.2
        ///<summary>Calculates the incomplete gamma function evaluated at A and X.</summary>
        ///<param name="a">the parameter of the gamma distribution.</param>
        ///<param name="x">the integration end point.</param>
        ///<returns>the incomplete gamma function evaluated at A and X.</returns>
        ///<remarks>The method returns <c>Double.NaN</c> if X &lt; 0 or A &lt;= 0.</remarks>
        ///<exception cref="ArithmeticException">if the the method is unable to converge to a value.</exception>
        public static double IncompleteGamma(double a, double x){
            if (Double.IsNaN(a) || Double.IsNaN(x) || (a <= 0.0) || (x < 0.0)) {
                return Double.NaN;
            } else if ( x==0 ) {
                return 0;
            } else if ( Double.IsPositiveInfinity(x) ) {
                return 1;
            } else {
                double an = 1.0/a;
                double sum = an;
                for (int n = 1; n <= MAX_ITERATIONS; ++n) {
                    an *= x / (a+n);
                    sum += an;
                    if (System.Math.Abs(an) < Double.Epsilon) {
                        return sum * System.Math.Exp(-x + a * System.Math.Log(x) - LnGamma(a));
                    }
                }
                throw new System.ArithmeticException("Unable to calculate Incomplete Gamma.");
            }
        }

        ///<summary>Calculates the beta function evaluated at P and Q.</summary>
        ///<param name="p">value to evaluate.</param>
        ///<param name="q">value to evaluate.</param>
        ///<returns>the beta function evaluated at P and Q.</returns>
        ///<remarks>The method returns <c>Double.NaN</c> if P,Q, P+Q &lt;= 0.</remarks>
        public static double Beta(double p, double q){
            return System.Math.Exp( LnGamma(p)+LnGamma(q)-LnGamma(p+q) );
        }

        ///<summary>Calculates the error function evaluated at X.</summary>
        ///<param name="x">value to evaluate.</param>
        ///<returns>the error function evaluated at X.</returns>
        public static double Error(double x){
            if ( Double.IsNaN(x) ) {
                return Double.NaN;
            } else if ( x == 0) {
                return 0;
            } else if ( Double.IsPositiveInfinity(x) ) {
                return 1;
            } else if ( Double.IsNegativeInfinity(x) ) {
                return -1;
            } else {
                return(x < 0) ? -IncompleteGamma(.5, x*x) : IncompleteGamma(.5, x*x); 
            }
        }

        ///<summary>Calculates the inverse error function evaluated at X.</summary>
        ///<param name="x">value to evaluate.</param>
        ///<returns>the error function evaluated at X.</returns>
        public static double InverseError(double x){
            if ( Double.IsNaN(x) ) {
                return Double.NaN;
            } else if ( x == 0) {
                return 0;
            } else if ( x == 1 ) {
                return Double.PositiveInfinity;
            } else {
                return(x < 0) ? -IncompleteGamma(.5, x*x) : 1-IncompleteGamma(.5, x*x); 
            }
        }

        ///<summary>Calculates the factorial of the integer N.</summary>
        ///<param name="n">value to evaluate.</param>
        ///<returns>the factorial of the integer N.</returns>
        public static double Factorial( int n ){
            if ( n < 0 ) {
                return Double.NaN;
            } else if ( n <= 20 ) {
                return FactorialTable[n];
            } else {
                return System.Math.Exp(LnGamma(n+1));
            }
        }

        ///<summary>Calculates the natural log of factorial of the integer N.</summary>
        ///<param name="n">value to evaluate.</param>
        ///<returns>the factorial of the integer N.</returns>
        ///<remarks>Returns 0 for N &lt;= 1. </remarks>
        public static double LnFactorial( int n ){
            if ( n < 0 ) {
                return Double.NaN;
            } else if ( n <= 1 ) {
                return 0.0;
            } else if ( n <= 20 ) {
                return System.Math.Log(FactorialTable[n]);
            } else {
                return LnGamma(n+1);
            }
        }

        ///<summary>Calculates the binomial function evaluated at N and K.</summary>
        ///<param name="n">value to evaluate.</param>
        ///<param name="k">value to evaluate.</param>
        ///<returns>the binomial function evaluated at N and K.</returns>
        ///<remarks>The method returns <c>Double.NaN</c> if N,K &lt; 0.</remarks>
        public static double Binomial( int n, int k ){
            if ( k == 0 ) {
                return 1;
            } else if ( k == 1) {
                return n;
            } else if ( k == n) {
                return 1;
            } else if ( k == (n+1) ) {
                return 0;
            } else if ( k == 2) {
                return(n-1)*n/2.0;
            } else if ( k == 3) {
                return(n-2)*(n-1)*n/6.0;
            } else if ( k == 4) {
                return(n-3)*(n-2)*(n-1)*n/24.0;
            } else {
                return System.Math.Floor(0.5+System.Math.Exp( LnFactorial(n)- LnFactorial(k)- LnFactorial(n-k)));
            }
        }
    }
}


